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ABSTRACT 

Motivated by a growing concern that masses of circumstellar disks may have been systematically 
underestimated by conventional observational methods, we present a numerical hydrodynamics study 
of time-averaged disk masses ((Afd)) around low-mass Class 0, Class I, and Class II objects. Mean 
disk masses (Md) are then calculated by weighting the time-averaged disk masses according to the 
corresponding stellar masses using a power-law weight function with a slope typical for the Kroupa 
initial mass function of stars. Two distinct types of disks are considered: self-gravitating disks, 
in which mass and angular momentum are redistributed exclusively by gravitational torques, and 
viscous disks, in which both the gravitational and viscous torques are at work. We find that self- 
gravitating disks have mean masses that are slowly increasing along the sequence of stellar evolution 
phases. More specifically. Class O/I/II self-gravitating disks have mean masses — 0.09, 0.10, and 
0.12 Mq, respectively. Viscous disks have similar mean masses (Md = 0.10 — 0.11 A/©) in the Class 
O/I phases but almost a factor of 2 lower mean mass in the Class II phase (Md,cii = 0.06 Mq). In 
each evolution phase, time-averaged disk masses show a large scatter around the mean value. Our 
obtained mean disk masses are larger than those recently derived by Andrews & Williams and Brown 
et al., regardless of the physical mechanisms of mass transport in the disk. The difference is especially 
large for Class II disks, for which we find Md,cii = 0.06 — 0.12 Mq but Andrews and Williams 
report median masses of order 3 x 10"'^ Mq. When Class O/I/II systems are considered altogether, a 
least-squares best fit yields the following relation between the time-averaged disk and stellar masses, 
(Md) = (0.2 ± 0.05) (M*)^-^='=°-^^. The dependence of (Md) on (M*) becomes progressively steeper 
along the sequence of stellar evolution phases, with exponents 0.7 ± 0.2, 1.3 ± 0.15, and 2.2 ± 0.2 for 
Class 0, Class I, and Class II systems, respectively. 

Subject headings: circumstellar matter — planetary systems: protoplanetary disks — hydrodynamics 
— ISM: clouds — stars: formation 



1. INTRODUCTION 

It has now become evident that disks of gas and dust 
are present from the earliest phases of stellar evolution 
(Class and Class I) and last for at least several million 
years into the late Class II phase. Disks are observed 
or inferred around most T Tauri stars and even around 
brown dwarfs. The evidence for disks around Class 
and Class I sources is more indirect. In this early phase 
of stellar evolution, the protostar/disk system is deeply 
embedded in an envelope - a remnant of the cloud core 
from which the protostar is forming. Nevert heless, recent 
observations bv'Andre ws fc Willi ams' ('2005) suggest that 
Class I disks have a larger median mass than that of Class 
II disks. 

In spite of a considerable progress in the detection 
of disks around young stellar objects (YSOs), an accu- 
rate determination of disk masses is still challenging. It 
is difficult to directly determine disk masses from the 
spectral lines of molecular species because the bright- 
est, easily detectable lines (i.e., the rotational transi- 
tions of CO) are optically thick and likely to be severely 
depleted. Therefore, disk masses are usually inferred 
from analyzing the spectral energy distribution of YSOs 
from the mid-infrared through submillimeter ba nds (e.g. 
[Andrews k Williamd [2005t iBrown et al] [2000h . Such 
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measurements of disk masses suffer from large uncertain- 
ties in the normaliz ation of dust opacitie s and gas-to-dust 
ratios, which led Ha rtmann et al.l (j2006f ) to conclude that 
T Tauri disk masses have been systematically underesti- 
mated by conventional analyses. 

Another complication arises from poorly known phys- 
ical processes in the disk. A usual assumption of op- 
tically thin circumstellar disks may significantly un- 
derestimate disk masses, particularly for objects with 
larger flux densities. However, a self-consistent treat- 
ment of a non-negligible optical depth requires a knowl- 
edge of the radial gas surface density profile in the 
disk (|Andrews fc Williamsl l2005f ). which may depend 
significantly on the dominant physical mechanism of 
mass and angular momen tum redistribution in the disk 
(jVorobvov fc Basull2008H) . 

Given large uncertainties in the measurements of disk 
masses, numerical simulations of self-consistent forma- 
tion and evolution of circumstellar disks can provide 
valuable information on disk masses in the early em- 
bedded and late phases of YSO evolution. It has been 
shown in the past that the saturation of spiral grav- 
itational instabilities at a finite amplitude in a self- 
gravitating, Toomre-unstable disk allows for the steady 
transport of mass an d momentum , which eventually lim- 
its disk masses fe.g I Adams et allll989t IShu et al.l(l990l : 
iLaughlin et all 119971 . \199^ In this paper, we perform 
numerical simulations of the long-term evolution of self- 
consistently formed circumstellar disks around low-mass 
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stars (0.2 Mq < A/, < 2.0 Mq). We consider both the 
self-gravitating disks, in which radial transport of mass 
and angular momentum is done exclusively via gravita- 
tional torques, and viscous disks, which feature gravita- 
tional torques as well as viscous ones. We seek to deter- 
mine numerically the disk masses in the Class 0, Class I, 
and Class II phases of stellar evolution. 

The paper is organized as follows. The numerical 
methods and initial parameters of cloud cores are given 
in § [21 Our obtained masses for self-gravitating and vis- 
cous disks are presented in § [3] and § S) We compare our 
numerical results with observations in § [S] The model 
and numerical caveats are discussed in § [51 The main 
results are summarized in § [71 

2. DESCRIPTION OF NUMERICAL METHODS 

We use the thin-disk approximation to compute the 
evolution of rotating, gravitationally bound cloud cores. 
This allows efficient calculation of the long-term evolu- 
tion of a large number of models. We start our numerical 
integration in the pre-stellar phase, which is character- 
ized by a collapsing starless cloud core, and continue into 
the main accretion phase, which sees the formation of a 
central star and circumstellar disk. We cover all major 
phases of the evolution of a YSO, starting from its for- 
mation and ending with the T Tauri phase. The integra- 
tion ends when the age of the central star is about three 
million years. In some models we extend the integra- 
tion up to 5 Myr. We emphasize that circumstellar disks 
are formed self-consistently in our numerical simulations, 
rather than being introduced as an initial parameter of 
the model. 

Once the disk is formed, its mass is determined by an 
interplay between the efficiency of the mass and angular 
momentum transport in the disk"^ and the infall rate of 
matter from the surrounding envelope onto the disk. At 
the time of disk formation, the infall rates take values 
between 1.2 x 10"^ Mq yr~^ and 7.0 x lO"*^ Mq yi'^ 
(measured at 600 AU) for the least and most massive 
cloud cores, respectively, but they show a fast decline 
with time. These values (and strong time variation) are 
consistent with the infall rates derived by Klessen (2001) 
using numerical models that follow molecular cloud evo- 
lution from turbulent fragmentation toward the forma- 
tion of stellar clusters. We note that once the disk is 
formed, the infall rate of matter from the envelope onto 
the disk is not necessarily the same as the mass accretion 
rate from the disk onto the protostar. While the for- 
mer shows a fast decline with time, the latter is usually 
characterized by a much slower decline and has a strong 
dependence on the stellar mass (IVorobvov fc Bas"ull2007l 
HoOSa). 

We use two numerical approaches: a basic approach 
that accounts for the radial transport due to gravita- 
tional toques and a viscous approach that accounts for 
the radial transport due to both the gravitational torques 
and viscosity. Gravitational torques are known to effi- 
ciently redistribute m ass and angular momentum in cir- 
cumstellar disks fe.g. iLodato fc Ricel[200l l2005f l. They 
were shown to play an important role in driving the 

^ In fact, disks may also transport angular momentum to the 
external environment due to magnetic braking. This effect will be 
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FU-Ori-lik e bursts in the early embedded phase of disk 
evolution (jVorobvov fc Basull2005bL l2006l ). In the late 
disk evolution, negative gravitational torques associated 
with low-amplitude azimuthal density perturbations in 
the disk can drive mass accretion rates that are consistent 
with those me a.sured in the interrnediate and u pper-mass 
T Tauri stars (jVorobvov fc Basull2007ll2008a[) . 

2.1. Basic numerical approach 

In the basic numerical approach, the collapse of a cloud 
core and subsequent evolution of a star/disk system is 
carried out by solving the basic equations of mass and 
mo mentum transport in the thin-disk approximation (see 
e.g. IVorobvov fc Basull2006D 



as 

'dt 
dt 



(1) 
(2) 



where S is the mass surface density, V = Pdz is 
the vertically integrated form of the gas pressure P, 
Z is the radially and azimuthally varying vertical scale 
height, Vp = VrV + v^cf) is the velocity in the disk plane, 
9p = drf" + is the gravitational acceleration in the 

disk plane, and Vp = rd/dr + (fir^^d/dcj) is the gradient 
along the planar coordinates of the disk. The gravita- 
tional accele ration g^, is found b y solving for the Poisson 
integral fsee IVorobvov fc Basul l20G6V The fact that we 
account for the disk self-gravity means that gravitational 
torques arise self-consistently in our numerical simula- 
tions and not imitated by some means of a-viscosity. 

Taking into account the complexity of gas thermody- 
namics in circumstellar disks (see §[6|), we have adopted a 
barotropic equation of state that closes equations ([T]) and 
^ and makes a transition from isothermal to adiabatic 
evolution at S = Scr = 36.2 g cm~^ 



r = cij: 



(3) 



where Cs is the isothermal sound speed, the value of 
which is set equal to that of the initial cloud core, and 
7 = 1.4. Equation ([3]), though neglecting detailed cool- 
ing and heating processes, was shown to reproduce to 
a first approximation the radial temperature gradients 
in the disk (Vorobyov fc Basu 2007) and the density- 
temperature relation fo r collapsing cloud cores derived by 
iMasunaga fc Inutsukal feOOO.) using a detailed radi ation 
hydrodynamics simulation (jVorobvov fc Basull2006D . 

The vertical scale height Z{r,(j),t) is determined as- 
suming the local hydrostatic equilibri um in the gravita- 
tional field of a disk and central star (jVorobvov fc Basu! 
i2008b) . The relevant formulas are given in the Appendix. 

2.2. Viscous numerical approach 

Viscosity is another important mechanism of angu- 
lar momentum and mass redistribution in astrophysical 
disks. Most analytical and numerical studies of viscous 
evolution of thin circumstellar d isks have employed the 
standa rd axisymmteric model of iLvnden-Bell fc Pringld 
(jl974[ ). in which the surface density of a Keplerian disk 
evolves with time according to the following diffusion 
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equation 



as 



3_9 

r dr 



(4) 



where v is the kinematic viscosity. 

In the present paper, we take a more fundamental ap- 
proach and describe the effect of (yet unspecified) viscos- 
ity in terms of the classic viscous stress tensor 



n = 2/1 f Vt; - i(V • v)e 



(5) 



where is a symmetrized velocity gradient tensor, e is 
the unit tensor, and /i is the dynamical viscosity. This 
approach allows for a self-consistent treatment of both, 
self-gravity and viscosity, within the same numerical for- 
malism. The resulting mass and momentum transport 
equations in the viscous numerical approach are 



9S 

'dt 
dv 



p 



dt 



-V„P + Sg„ + (V-n) 



(6) 
(7) 



where V • 11 is the divergence of the rank-two viscous 
stress tensor 11. The relevant components of (V • 11)^ 
are given in the Appendix. Equations ([6]) and ([7]) are 
closed with the barotropic equation of state ^ . We note 
that the viscous approach accounts self-consistently for 
both the gravitational and viscous torques which may 
arise during numerical simulations. 

It is evident that the practical application of equa- 
tion ([7]) requires a knowledge of the dynamical viscos- 
ity fi of the disk. Unfortunately, our understanding of 
viscous processes in circumstellar disks is still incom- 
plete. We know that molecular (collisional) viscosity 
is most certainly too low to be of practical interest. 
Turbulence driven by the magneto-rotational instability 
(MRI) is a most promising source of viscosity at present 
(|Balbus fc Hawlev 1991), though other mechanisms can- 
not be ruled out completely. 

In this paper, we make no specific assumptions as 
to the source of viscosity and define the coefficient of 
dynamical viscosity using the usual a-prescription of 
IShakura fc SunvaevI (I1973D 

/i = a S Cs (8) 

where spatially and temporally uniform a is set equal to 
0.01 and Cg — \J dV /dY, is the effective sound speed. Our 
numerical simulations of embedded and T Tauri disks in- 
dicate that a ~ 0.001 — 0.01 yields disk sizes and radial 
slopes that are in genera l agrement with observations 
(|Vorobvov fc Basul[2008bf ). Lower values of a (< 10"^) 
have little effect on the evolution of self-gravitating disks, 
whereas substant ially higher valu e s (> 0.1) quickly de- 
stroy the disks. iHartmann et al.l (|1998f ) also predicted 
similar values for a by analyzing accretion rates in T 
Tauri disks. Since viscosity in our model is assumed to 
arise due to some physical processes in the disk, we keep 
a equal to zero during the early "pre-disk" phase of evo- 
lution and set a equal to 0.01 only when a circumstellar 
disk forms around a central star. 

2.3. Initial condition 



TABLE 1 

Parameters of models with rji = 1.2 x 10"'^ 



Model 


ro 


So 


Ho 


rout 




1 


1209 


0.13 


1.14 


7186 


0.7 


2 


1382 


0.12 


1.0 


8213 


0.8 


3 


1728 


0.093 


0.8 


10266 


0.98 


4 


2074 


0.077 


0.67 


12320 


1.18 


5 


2937 


0.055 


0.47 


17452 


1.67 


6 


4147 


0.039 


0.33 


24640 


2.36 


7 


5184 


0.031 


0.27 


30800 


2.95 



Note. — All distances are in AU, angular velocities in 
km s"-"^ pc""*^, surface densities in g cm~^, and masses in Mq. 



The initial radial distributions of surface density S and 
angular velocity in our model cloud cores are those 
characteristic of a collapsing axisymmctric magnetically 
supercritical core (|Basuill997D 




(9) 



(10) 



where ro is the radial scale length defined as Tq = 
kcl/{GT,o) and k = V^/n. These initial profiles are 
characterized by the important dimensionless free pa- 
rameter rj = ripTg/Cg and have the property that the 
asymptotic (r ^ rp) ratio of centrifugal to gravitational 
acceleration has mag nitude V^tj (see iBasul [19971) . The 
centrifugal radius of a mass shell initially located at ra- 
dius r is estimated to be rcf = p /{Gm) = \/2rir, where 
j — flr^ is the specific angular momentum. Since the en- 
closed mass m is a linear function of r at large radii, this 
also means that r^f oc m. The gas has a mean molecu- 
lar mass 2.33 mn and cloud cores are initially isothermal 
with temperature T = 10 K. 

We present results from three sets of models, each 
with a different value of rj. The standard model has 
T] ^ i]i ~ 1.2 X 10"'^ based on typical values Cg = 0.19 
km s^^. So = 0.12 g cm~^, and fig = 1.0 km s~^ pc~^. 
The outer radius is taken to be Tout = 0.04 pc, and the 
total cloud mass is 0.8 Mq. Other models with 77 = 771 
but different mass (outer radius) are generated by vary- 
ing ro and fio so that their product is constant. Note 
that, when ro is varied. So has to be changed accord- 
ingly. All clouds are characterized by the same ratio 
Tout I To ~ 6.0. To generate the second set of models, 
V = V2 = 2.3 X 10"'^, we set ^Iq = 1.4 km pc~^ and all 
other quantities the same as in the standard model with 
77 = ?7i . Models of varying mass are then generated in the 
same manner as for the 771 models. The third set of mod- 
els, with rj = ri3 = 3.4 x 10~^, are also obtained in this 
way, by first using ^lo — 1.7 km pc^^. Overall, there 
are 7 models with rj = iji, 13 models with 77 = 772 — 2771, 
and 12 with 77 = 773 ~ 3771. The range of initial cloud 
masses (Md) amongst our models is 0.3 Mq — 2.95 Mq. 
The parameters of our models are listed in Table [T] (771), 
Tabled (772), and Tabled (773). We note that our model 
values of ilo = (1.0— 1.7) km s~^ pc~^ are within a typi- 
cal rang e of velocity g r adien ts measured in dense starless 
cores bv lCaseUi etall (I200I . 
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TABLE 2 

Parameters of models with ri2 = 2.3 x 10"'' 



Model 


ro 


So 


Ho 


rout 


A/ci 


8 


622 


0.26 


3.1 


3696 


0.35 


9 


691 


0.23 


2.8 


4106 


0.4 


10 


864 


0.19 


2.24 


5133 


0.5 


11 


1037 


0.16 


1.87 


6160 


0.6 


12 


1210 


0.13 


1.6 


7187 


0.7 


13 


1417 


0.11 


1.4 


8213 


0.8 


14 


1728 


0.093 


1.12 


10267 


0.98 


15 


2454 


0.065 


0.79 


14579 


1.4 


16 


2765 


0.058 


0.7 


16428 


1.57 


17 


3456 


0.046 


0.56 


20533 


1.97 


18 


3802 


0.042 


0.51 


22587 


2.16 


19 


4147 


0.038 


0.47 


24640 


2.36 


20 


4838 


0.033 


0.4 


28747 


2.75 



Note. — All distances are in AU, angular velocities in 
km pc~^, surface densities in g cm~^, and masses in Mq. 



TABLE 3 

Parameters of models with ri2 = 3.4 x 10~^ 



Model 


ro 


So 


Uo 


Tout 


A/ci 


21 


518 


0.31 


4.53 


3080 


0.3 


22 


691 


0.23 


3.4 


4106 


0.4 


23 


864 


0.19 


2.7 


5133 


0.5 


24 


1037 


0.16 


2.26 


6160 


0.6 


25 


1210 


0.13 


1.94 


7187 


0.7 


26 


1417 


0.11 


1.7 


8213 


0.8 


27 


2073 


0.077 


1.13 


12320 


1.18 


28 


2420 


0.066 


0.97 


14373 


1.38 


29 


3283 


0.049 


0.72 


19506 


1.87 


30 


3802 


0.042 


0.62 


22587 


2.16 


31 


4147 


0.039 


0.57 


24640 


2.36 


32 


4493 


0.036 


0.52 


26693 


2.56 



Note. — All distances are in AU, angular velocities in 
km pc~^, surface densities in g cm~^, and masses in Mq. 



2.4. Numerical technique 

Hydrodynamic equations of the basic and viscous nu- 
merical models are solved in polar coordinates (r, 0) on a 
numerical grid with 128 x 128 points. We use the method 
of finite differences with a time-explicit, operator-split 
solution procedure. Advection is performed using the 
second-order van Leer scheme. The radial points are log- 
arithmically spaced. The innermost grid point is located 
at r = 5 AU, and the size of the first adjacent cell lies 
in the range between 0.26 AU (model 21) and 0.36 AU 
(model 7). It means that the ratio Ar/r of the cell size 
Ar to radius r is constant for a given cloud core and 
varies from 0.05 (model 21) to 0.07 (model 7). 

We introduce a "sink cell" at r < 5 AU, which repre- 
sents the central star plus some circumstellar disk ma- 
terial, and impose a free inflow inner boundary condi- 
tion. About 95 per cent of the material crossing the inner 
boundary lands into the star, the rest constitutes an in- 
ner circumstellar region of constant surface density. The 
dynamics of the inner region (r < 5 AU) is not computed 
but it contributes to the total gravitational potential of 
the system. We do not account for a possible mass loss 
from the sink cell due to stellar jets, since our model is 
two-dimensional and it is not clear how the mass ejection 
efficiency varies with stellar age. The outer boundary is 



such that the cloud has a constant mass and volume. 
More detail s on numerical technique s and relevant tests 
are given in IVorobvov fc Basul (|2006D . 

3. MASSES OF SELF-GRAVITATING DISKS 

In this section we present disk masses obtained using 
the basic numerical approach described in § 12.11 In the 
framework of this numerical model, disk masses are con- 
trolled by the rate of mass infall from the surrounding en- 
velope and the radial transport of mass and angular mo- 
mentum due to gravitational torques. No viscous torques 
are present in this case. 

An accurate determination of disk masses in numeri- 
cal simulations of collapsing cloud cores is not a trivial 
task. Self-consistently formed circumstellar disks have 
a wide range of masses and sizes, which are not known 
a priori. However, numerical and observational stud- 
ies of circumstellar disks indicate that the disk surface 
density is a declining function of radius. Therefore, we 
distinguish between disks and infalling envelopes using 
a critical surface density for the disk-to-envelope transi- 
tion, for which we choose a value of Str = 0.1 g cm~^. 
This choice is dictated by the fact that densest starless 
cores have surface densities only slightly lower than the 
adopted value of Etr- In addition, our numerical simula- 
tions indicate that self-gravitating disks have sharp outer 
edges and the gas densities of order 0.01 — 0.1 g cm~^ 
characterize a typical transition region between the disk 
and envelope (Vorobyov & Basu 2007j). 

To compare disk masses in three distinct phases of stel- 
lar evolution we need an evolutionary indicator to dis- 
tinguish between Class 0, Class I, and Class II phases. 
We use a classification of lAndre et alj (|1993D . who sug- 
gest that the transition between Class and Class I ob- 
jects occurs when about 50% of the initial cloud core is 
accreted onto the protostar/disk system. The Class II 
phase is consequently defined by the time when the in- 
falling envelope clears and its total mass drops below 10% 
of the initial cloud core mass M^i. It should be mentioned 
here that there exists no unique classification schem e for 
protostars. For instance. IVorobvov fc Basiil (|2005af l have 
proposed a classification scheme that hinges on a tempo- 
ral behaviour of bolometric luminosity Lboi- They iden- 
tify Class phase with a period of temporally increasing 
Lboi and Class I phase with a later period of decreasing 
Lboi- The peak in Lboi corresponds to the evolutionary 
time when 50 ± 15 per cent of the cloud core mass has 
been accreted by the protostar. Observers prefer to clas- 
sify proto stars by thei r spec tral energy distributio n s. For 
instance, iLadal |T98l and [Andrews &: WilUamsl (pOOSh 
use the values of power-law index n (defined by vFi, oc 
I/"-, where F^, is the infrared flux density at frequency 
v) to distinguish b etween Class I and Class II objects. 
lAndre et al.l ()1993f ) use the ratio of submillimeter to bolo- 
metric luminosity isubmm/^boi for the same purposes. 
Although these classifications are physically related, it is 
possible that they may differ from our adopted classifi- 
cation, especially if some of the gas is removed from the 
cloud core and does not accrete. We acknowledge that 
the difference in the existing classification schemes may 
systematically shift our results. 

3.1. Temporal evolution of self-gravitating disks 
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Fig. 1. — Disk masses (thick solid lines), stellar masses (thin solid 
lines), and envelope masses (dashed lines) obt ained in the frame- 
work of the basic numerical approach (see ii l2.1l l for model 8 (upper 
left), model 11 (upper right), model 16 (lower left), and model 20 
(lower right). All masses are calculated relative to the correspond- 
ing initial cloud core mass M^i. The horizontal axis shows time 
elapsed since the formation of a central star. The left/right verti- 
cal dotted lines mark the onset of Class I/II phases, respectively. 



We start by comparing the long-term evolution of disk 
masses in four sample models, which are chosen to repre- 
sent a full spectrum of the initial cloud core masses with 
a constant ratio rj = 2.3 x lO^'^ of centrifugal to grav- 
itational acceleration. Figure [1] shows the disk masses 
(thick solid lines), stellar masses (thin solid lines), and 
envelope masses (dashed lines) in model 8 (upper left, 
Mci = 0.35 Mq), model 11 (upper right, Md = 0.6 Mq), 
model 16 (lower left, M^i = 1.57 Mq), and model 20 
(lower right, M^i — 2.75 Mq). The horizontal axis is the 
time elapsed since the formation of a central star. All 
masses are calculated relative to the corresponding ini- 
tial cloud core mass Md- The vertical dotted lines mark 
the onset of Class I (left) and Class II (right) phases. 

Our numerical simulations demonstrate that cloud 
cores with constant rj (but different size) form disks 
roughly at the same physical time after the formation 
of a central star but in distinct stellar evolution phases. 
For instance, model 8 starts to build a disk in the late 
Class I phase, while model 20 does that in the midst of 
Class phase. Cloud cores of greater mass tend to form 
disks in the earlier phase of stellar evolution than their 
low-mass counterparts. 

The upper panel of Fig. [T] identifies cases when Class 
stars have no disks associated with them. We emphasize 
here that due to the use of the sink cell in our numerical 
code we can resolve only those disks whose outer radii are 
larger than 5 AU. It means that even though circumstel- 
lar disks do not form around some Class objects in our 
numerical simulations, as in models 8 and 11 (top panels 
in Fig. [3]), one may suppose that such disks still exist but 
their size is simply smaller than that of the sink cell. Our 
test runs with a sink cell set to 0.5 AU (instead of 5 AU) 
confirm that disks indeed forms earlier in the evolution 
but quickly expands to 5 AU and beyond. Hence, the 
absence of disks around some Class objects may be to 
some extent caused by a finite-size sink cell. However, it 
is still possible that some Class objects have no disks 
associated with them. This is particularly true for those 
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Fig. 2. — Time-averaged disk masses (right column) and 
time-averaged disk-to-star mass ratios (left column) versus time- 
averaged stellar masses. Open triangles, filled squares, and plus 
signs show the data for Class 0, Class I, and Class II systems, re- 
spectively. In particular, top panels show the data for models with 
•qi = 1.2 X 10~^, whereas middle and bottom panels show the data 
for models with ri2 ~ 2r)i and rj^ ~ 3r)i , respectively. 

objects that develop from cloud cores with low rates of 
rotation, since in this case a (substantial) portion of the 
disk material will be characterized by the centrifugal ra- 
dius (estimated as = p /{Gm)) that is smaller than 
the radius of the protostar, 4 Rq. A numerical code with 
realistic treatment of the protostar formation is needed 
to accurately address this issue. 

Another interesting feature of self-gravitating disks 
that can be seen in Fig. [1] is that the disk-to-star mass 
ratio never exceeds some characteristic value, approx- 
imately 0.35 — 0.4, irrespective of the initial cloud core 
mass. This is rather counterintuitive. For instance, mod- 
els 16 and 20 (bottom panels in Fig. [1]) form disks in 
the Class phase, which is characterized by envelope 
masses that are considerably greater than those of the 
protostars. As a consequence, one may expect the for- 
mation of a disk with mass that is at least comparable 
to or even greater than that of a protostar. It turns 
out, however, that circumstellar disks that form in the 
Class or early Class I phase around stars with mass 
M* > 0.6 Mq develop vigourous gravitational instability 
- a very efficient means of inward mass transport that 
helps keep the disk mass well below that of the proto- 
star. Sharp drops in the disk mass (or equivalent surges 
in the stellar mass) seen in the upper right and bottom 
panels o f Figure [1] are a manifestation o f this process 
seee.g. fv^obvov fc Basull2005U [20061 ). iKratter et all 
2008( ) also predict that disks around stars with mass 
greater than 1.0 Mq are expected to be vigorously grav- 
itationally unstable in the early embedded phase of disk 
evolution. The late evolution phase (Class II) sees only 
an insignificant decline of disk mass with time. 
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3.2. Disk masses versus stellar masses 

We use 32 models, the parameters of which are listed 
in Tables [l]|3l to analyse the statistical relations between 
time-averaged disk and stellar masses in the Class O/I/II 
evolution phases. Time averaging is performed sepa- 
rately for each evolution phase over the duration of the 
phase. The age of the oldest disk in our sample is about 
3 Myr. Since disk masses in the Class II phase have a 
tendency to decline with time and the actual lifetime of 
Class II objects may be longer, we expect that we might 
have somewhat overestimated disk masses around Class 
II objects. This, obviously, does not effect our estimates 
of Class O/I disk masses. 

Figure [5] shows time-averaged disk masses (right col- 
umn) and time-averaged disk-to-star mass ratios (left 
column) versus time-averaged stellar masses for Class 
(open triangles). Class I (filled squares) and Class II ob- 
jects (plus signs). Several interesting conclusions can be 
drawn by analysing the figure. 

1. Time-averaged disk and stellar masses ((Afd) and 
(M,), respectively) in models with the same ratio 
of rotational to gravitational acceleration fall onto 
a unique evolutionary track in the (Md) — (M,) 
diagram. 

2. Class objects occupy the lower-left part of each 
evolutionary track. No stars with (M*) > 0.9 Mq 
have Class disks (but they have Class I/II disks). 
In other words, only low-mass stars (within our 
range of interest, 0.2 Mq < (M*) < 2.0 Mq) can 
harbour Class disks. 

3. Stars of equal mass have disks with similar masses, 
regardless of the stellar evolution phases. 

4. In each stellar evolution phase, disk-to-star mass 
ratios tend to have greater values for stars of 
greater mass. However, there is a clear saturation 
effect for stars with masses greater than 1.0 Mq 
- disk masses never grow above 40% of the stellar 
masses, even for models with the largest values of 
rj = r]3 — 3.4 x 10~^. The saturation of disk-to- 
star mass ratios is caused by the onset of vigorous 
gravitational instability in circumstellar disks that 
form in the Class and early Class I phase. 

We summarize the main properties of our model self- 
gravitating disks in Table ID In order to facilitate the 
comparison of our numerical results with observations, 
we calculate mean disk masses (in each evolution phase) 
by weighting the time-averaged disk masses according to 
the corresponding stellar masses. The relative impor- 
tance of the stellar masses is calculated using the follow- 
ing power-law weight function 

^Ki(Ai*;j-| ^^^,^^^_2.3 if (Af^) > 0.5 Mo . ^^^^ 

The slope of JFk ((-^*)) is typ ical for the K roupa initial 
mass function (IMF) of stars (|Kroupall200l . The mean 
disk masses are then calculated as follows 



Time {Myr) 



1 0.01 




Md,ph 



Eti(^lph)-^K((M;p,)) 

Eti^K((M:,p,; 



(12) 



Fig. 3. — The same as Figure [T] but obtained in th e framework 
of the viscous numerical approach described in § 12.21 

where {M^ and (M* pj^) are the time-averaged disk 
and stellar masses in the i-th model (see Fig. [2]) and in- 
dex "ph" stands for Class (CO), Class I (CI), or Class 
II (CII) stellar evolution phases. The summation is per- 
formed separately for each evolution phase. In partic- 
ular, there are TV = 20 models that have Class disks 
and = 31 models that have Class I disks. All 32 
models have Class II disks, of course. The ratio of the 
normalization constants A/B ~ 2 needed to evaluate 
equation (jl2p is found using the continuity condition at 
(M,) = 0.5 Mq. 

Table 2] indicates that the mean masses slightly in- 
crease along the stellar evolution sequence from Md,co = 
0.09 Mq around Class objects to Md,cii = 0.12 Mq 
around Class II objects. The minimum and maximum 
time-averaged disk masses ((M™p'J^) and {M™^f^), respec- 
tively) in each evolution phase are shown in the last three 
columns of Table ID It is seen that both the Class I and 
Class II disks have a similar range of masses, while Class 
disks have a somewhat narrower range. The smallest 
time-averaged mass of a self-gravitating disk found in our 
numerical simulations is (Md) — 0.02 Mq. 

4. MASSES OF VISCOUS DISKS 

In this section we present disk masses obtained in the 
framework of the viscous numerical approach described 
in § 12.21 We seek to determine the effect of viscosity, 
and associated viscous torques, on the masses of self- 
consistently formed circumstellar disks. Our comparison 
study of viscous versus self-gravitating disks have shown 
that viscous disks may lack a sharp transition bound- 
ary be tween the disk and the envelope (|Vorobvov &: Basul 
l2008bD . This smearing of the disk's outer physical bound- 
ary is particularly pronounced in disks of low mass and in 
the late Class II phase. In this situation, it is somewhat 
arbitrary to adopt a value of Str = 0.1 g cm~^ for the 
disk-to-envelope transition. Nevertheless, we have recal- 
culated disk masses with Etr = 0.01 g cm~^ and found 
that this can increase masses of Class II disks by 25% at 
most. Class O/I disks are unaffected by the value of Etr, 
since they usually have sharp outer physical boundaries. 

4.1. Temporal evolution of viscous disks 



Disk masses around young stellar objects 



7 



TABLE 4 
Summary of disk properties 



disk type A/d,co 


A/d,ci 


A/d.cii 




(A^d^co)) ((A-^dTi) 






(A^d":cn)) 


self-gravitating 0.09 
viscous 0.10 


0.10 
0.11 


0.12 
0.06 


(0.04 
(0.05 


: 0.2) (0.02 
0.19) (0.02 


: 0.5) 
0.55) 


(0.03 
(0.004 


0.7) 
: 0.6) 



Note. — Mean masses Afd.cOi Afd.cii ^^^d Afd.cil °f Class 0, Class I, and Class II disks, respectively, are calculated according to 
equation II12I I. Angle brackets () denote time averaging over the duration of a particular stellar evolution phase. Indices "min" and "max" 
refer to minimum and maximum time-averaged disk masses, respectively. All disk masses are in Mq. 



In order to investigate the temporal evolution of vis- 
cous disks, we consider the same four sample model cloud 
cores as in § 13.11 Figure [3] presents disk masses (thick 
solid lines), stellar masses (thin solid lines), and envelope 
masses (dashed lines) in model 8 (upper left), model 11 
(upper right), model 16 (lower left), and model 20 (lower 
right). The horizontal axis is time elapsed since the for- 
mation of a central star. All masses are calculated rel- 
ative to the corresponding initial cloud core mass Md- 
The vertical dotted lines mark the onset of Class I phase 
(left line) and Class II phase (right line). 

The comparison of Figures [T] and [3] reveals a major dif- 
ference between self-gravitating and viscous disks - the 
latter have considerably smaller masses in the late Class 
II phase than the former. For instance, model 8 in Fig- 
ure[l]has a final {t = 5 Myr) stellar mass = 0.33 Mq 
and disk mass Md — 2.5 x 10^^ Mp, while the same 
model in Figure[3]has Md — 8.5 x 10~ Mq. Such a large 
contrast, however, diminishes for stars of greater mass. 
For instance, model 16 in Figure [T] has M* = 1.57 Mq 
and Md = 0.37 Mq, while the same model in Figure [3] 
has Md = 0.11 Mq. On the other hand, both the vis- 
cous and self-gravitating disks have comparable masses 
in the Class O/I phases. Even the early Class II phase 
sees little difference between the masses of viscous and 
self-gravitating disks, provided that the disks are formed 
from molecular cloud cores of similar mass. Figure [3] in- 
dicates that a small portion of the viscous disk material 
returns to the envelope in the late disk evolution. This 
is because we have imposed a fixed gas column density 
threshold for the disk-to-envelope transition, but viscous 
disks expand in the late evolution due to the action of vis- 
cous torq ues, whic h are in fact positive in the outer disk 
regions (Vorojjyw & Basu 2008b) . This numerical effect 
may slightly reduce the resulted viscous disk masses in 
the Class II phase. 

Why does the late stellar evolution phase reveal such 
a noticeable difference in the disk masses (and, conse- 
quently, in disk-to-star mass ratios) between viscous and 
self-gravitating disks and why is there little difference in 
the earlier phases? The reason for that can be under- 
stood if we consider the temporal evolution of viscous 
and gravitational torques in the disk. Figure [4] shows 
total negative torques due to disk self-gravity (7^, solid 
lines) and viscosity (7^, dashed lines) in our four typical 
models: model 8 (upper left), model 11 (upper right), 
model 16 (lower left), and model 20 (lower right). The 
horizontal axis shows time elapsed since the formation of 
a central star. The total negative torque is computed by 
summing up local negative torques r(r, 0) in each com- 
putational zone occupied by the disk. In particular, 7^ is 
computed as a sum of rg(r, 0) = ~m{r, (/)) d^{r, (j)) /d(f), 
where m(r, 0) and <I>(r, (j)) are the gas mass and gravi- 
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Fig. 4. — Total negative gravitational torques (solid lines) and to- 
tal negative viscous torques (dashed lines) versus time elapsed since 
the formation of a central star in model 8 (upper left), model 11 
(upper right), model 16 (lower left), and model 20 (lower right). 
The vertical dotted lines mark the onset of Class I phase (left line) 
and Class II phase (right line). The torques are in units of 8.7x 10^" 
g cm s . 

tational potential in a cell with polar coordinates (r, 0), 
respectively. The total negative viscous torque is calcu- 
lated in a similar manner by summing up all local neg- 
ative viscous torques T^{r, 4>) — r {y ■ n)^ S{r, 4>), where 
S{r, (/)) is the surface area occupied by a cell with polar 
coordinates (r, (j)). The local torques Tg{r, cj)) and T-y{r, cj)) 
are actually the local azimuthal components of the cor- 
responding forces, acting on a fluid element with mass 
m(r, 0), multiplied by the arm r. 

Figure d] indicates that gravitational torques always 
prevail over viscous torques in the Class and Class I 
phases. As a result, the disk mass in these phases is de- 
termined by the interplay between the mass load from 
an infalling envelope and the rate of mass and angular 
momentum transport in the disk due to gravitational 
torques. In the Class II phase, the strength of gravi- 
tational torques diminishes, partly due to a stabilizing 
influence of a growing central star and partly due to the 
exhausted mass reservoir in the infalling envelope. As a 
result, viscous torques, the strength of which shows only 
a mild decline with time, take a leading role in the Class 
II phase. They act to further decrease the disk mass 
when gravitational torques fail to do so. 

The actual time in the Class II phase when %, become 
greater than 7^ is distinct for models with different ini- 
tial cloud core masses (but similar ratios of rotational 
to gravitational acceleration rf). Low-mass cloud cores 
form disks late in the evolution (see e.g. model 8 in 
Fig. [3]). The resulting disks have low masses, merely due 
to the fact that most of the cloud core material has al- 
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Fig. 5. — The same as Figure [2] only obtained in th e framework 
of the viscous numerical approach described in ii l2.2l 

ready been accreted directly onto the protostar in the 
pre-disk phase. As a result, gravitational instability is 
underdeveloped and viscous torques quickly take over the 
gravitational ones in the beginning of the Class II phase. 
Conversely, high-mass cloud cores form disks early in the 
evolution (see e.g. model 16 and 20 in Fig. [3]). The re- 
sulting disks are characterized by considerable masses. 
As a consequence, gravitational torques prevail over vis- 
cous torques throughout a considerable portion of the 
Class II phase. The in terested reader is referred to 
IVorobvov fc Basul (|2008al lbl) for a detailed discussion on 
gravitational torques and appearances of self-gravitating 
and viscous disks. 

4.2. Relation between disk and stellar masses 

In this section we derive statistical relations between 
time-averaged disk masses ((Md)), time-averaged stellar 
masses ((M*)) , and time-averaged disk-to-star mass ra- 
tios ((0) for viscous disks in each stellar evolution phase. 
Time averaging is performed separately in each evolu- 
tion phase over the duration of the phase. Figure [5] 
shows (Afd) (right column) and (^) (left column) ver- 
sus (M*) for Class objects (open triangles). Class I 
objects (filled squares) and Class II objects (plus signs). 
The top-to-bottom rows show the data for models with 
r]i = 1.2 X 10~^ T]2 = 2.3 X 10"^, and 773 = 3.4 x 10~^ 
respectively. 

The comparison of Figures [5] and [H] reveals a distinct 
feature of viscous disks around objects of equal stellar 
mass - Class II disk have systematically lower masses 
than their Class O/I counterparts. This is in sharp con- 
trast to self-gravitating disks studied in § [3] where we 
have found that Class O/I/II objects with equal stellar 
mass possess self-gravitating disks of essentially similar 
mass (see Fig. [5]). The difference in masses between Class 
II and Class O/I viscous disks is considerable for low- 



mass stars (up to a factor of 10), but it diminishes for 
stars of greater mass. However, we have to keep in mind 
that low-mass stars are statistically more important than 
their high-mass counterparts. 

The dependence of {£,) on (Af,) in Figure [5] is similar 
to that of Figure [2] - stars of greater mass tend to have 
greater disk-to-star mass ratios. However, as for the case 
of disk masses, Class II viscous disks have disk-to-star 
mass ratios that are lower than those of self-gravitating 
disks. For instance, the smallest (^) in self-gravitation 
disks (Fig. [2]) is 0.05, while viscous disks have (^) as 
small as 0.01. 

It is in teresting to co r npare our results with the recent 
work of iKratter et al.l (j2008f ). who use semi-analytical 
models to describe the time evolution of embedded, ac- 
creting protostellar disks. As in our model, they ac- 
count for disk self-gravity and MRI-induced viscosity but 
also include the effect of stellar and envelope irradiation. 
They obtain maximum disk-to-star mass ratios in the 
Class I phases of order 30% — 40%, in excellent agree- 
ment with our results. The disk mass in their fiducial 
1 Mq model shows a temporal behaviour similar to that 
shown in the lower-left panel of Fig. [3]- it grows in the 
Class phase, reaches a maximum in Class I phase, and 
declines afterwards. They also predict that higher mass 
stars have relatively more massive disks and this trend is 
expected to extend to stellar masses much greater than 
those considered in our work (< 2.0 Mq). 

We calculate the mean masses of viscous disks us- 
ing equation (|12p and summarize the main properties 
of our viscous disks in Table [H It is evident that the 
range of time-averaged masses and the values of mean 
masses are quite similar in both the self-gravitating 
and viscous disks. The latter tend to have slightly 
larger mean disk masses in the Class O/I phases, which 
is explained by the fact that viscous torques tend to 
oppose gravitational one s in the early disk evolution 
(jVorobvov fc Basul l2008bD . The only significant differ- 
ence comes from Class II viscous disks, which have a 
factor of two lower mean mass. Men = 0.06 Mq, than 
that of self-gravitating disks. Class II viscous disks also 
feature the lowest time-averaged disk mass found in our 
numerical simulations, (M^j^^-^) = 0.004 Mq. 

Finally, we seek to determine the relation between 
time-averaged disk and stellar masses in each stellar evo- 
lution phase separately and for the three phases taken 
altogether. Figure [6] shows (Md) versus (M*) for our 
32 model cloud cores. In particular, the open triangles, 
filled squares, and plus signs show the data for the Class 
0, Class I, and Class II phases, respectively. It is evident 
that the upper-mass stars have a considerably narrower 
scatter in disk masses than the lower-mass stars. The 
least-squares best fit to all data in Figure [6] (solid fine) 
yields the following relation between time-averaged stel- 
lar and disk masses (in Mq) for the three phases taken 
altogether 

(Md) = (0.2 ± 0.05) (M,)i-3±o-i5. (13) 

The relation between time-averaged disk and stellar 
masses in each stellar evolution phase is, however, quite 
different from that expressed by equation (fT5|) . The least- 
squares best fits performed separately for each evolution 
phase yield the following relations 

(Md,co> = (0.2 ± 0.05) (M,,co>°-'±°-^ (14) 
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Fig. 6. — Time-averaged disk masses versus time-a vera ged stel- 
lar masses for 32 models of our sample (see Tables rH3l l. Open 
triangles, filled squares, and plus signs show the data for Class 0, 
Class I, and Class II systems, respectively. The solid line gives a 
least-squares best fit to all data in the plot. 

(Md,ci) = (0.3 ± 0.04) (M.,cI)^■^±°■^^ (15) 

(Md,cii) = (0.2 ± 0.08) (M,,cii>'-'^" ', (16) 

where indices CO, CI, and CII refer to the Class 0, 
Class I, and Class II phases, respectively. The depen- 
dence of disk masses on stellar masses becomes pro- 
gressively steeper along the sequence of stellar evolution 
phases. 

The above relations between time-averaged stellar and 
disk masses were derived for viscous disks, in which 
both viscosity and self-gravity are at work. Purely self- 
gravitating disks (zero viscosity) have a similar rela- 
tion between time-averaged disk and stellar masses when 
all stellar evolution phases are considered altogether, 
namely equation (|13p. However, when each phase is 
taken separately, the exponents are different from those 
derived for viscous disks in equations p^ - pB]) . To avoid 
confusion, we provide here the relations for viscous disks 
only, since we believe that some sort of viscosity should 
be present in circumstellar disks. 

Observations provide conflicting evidence as to the re- 
lation between disk and st ellar masses in YSOs. For in- 
stance, iNatta et al.l (|200lD found a marginal correlation 
between disk and stellar masses, a lbeit with a substan- 
tial di spersion. On the other hand, 'Andrews, & Williams 
(|2005f ) claimed no correlation. The correlation in Fig- 
ure[6]can be broken if a large population of objects occu- 
pies both the upper-left and lower-right portions of the 
(Aid) — {M^,) diagram. While the latter is feasible due 
to the incompleteness of our sample of model cloud cores 
(we have not considered sufficiently low values of the ra- 
tio Tj) , the former is unlikely due to the saturation of disk 
masses discus sed in the context o f Fig ure [H and earlier 
considered byEdamseFaD (|1989( l and lShu et all (|1990l ) 
. No stars can harbour disks with masses comparable to 
or greater than the stellar mass! It is likely that the un- 
certainty in disk measurements may introduce a consid- 
erable scatter to the observational data which precludes 
the correlation analysis. 

5. DISCUSSION 

There exist a growing concern that disk masses around 
YSOs may be systematically underestimated by conven- 



tional observational methods. The simplest argument in 
favour of this conjecture is a high frequency of detection 
of Jupiter-like exosolar planets. Both the core accretion 
and gravitational instability models require disk masses 
of order 0.1 Mq for Jupiter- like planets to form, which is 
about an order of magnitude higher tha n medi an disk 
masses inferred by [Andrews fc Williams! (|2005[ ) (here- 
after, AW). Another argument in favour of fairly mas- 
sive disks is the fact that resolved disks show rich non- 
axisymmetric structures, such as multi ple spiral arms or 
arcs in the case o f AB Aurigae fiFukagawa et al. 120041 : 
iCradvet al.l[l999l ) and HD 100546 (iGradv et al.lbooJ) . 
Such structures are difficult to sustain in low-mass disks. 
To remind the reader, a companion star or giant planet 
would most likely trigger a two-armed spiral response in 
the disk of a target and not a flocculent multi-arm struc- 
ture. 

As was discussed in the Introduction, measurements of 
disk masses suffer from the uncertainties in dust opacities 
and disk radial structure, which could l ead to a substan- 
tial u nderestimate of disk masses (e.g. iHartmann et al.l 
l2006f ). Alternatively, the disk may mask itself, especially 
late in the evolution, by locking a significant amount of 
its mass in the form of large dust grains or planetesi- 
mals which are inefficient emitters at submillimeter wave- 
lengths (|And rews & WiUiams 2005). 

Our numerica l mode ling supports the supposition of 
IHartmann et al.l (120061 ) that circumstellar disks are more 
massive than currently inferred from observations. The 
mean masses of disks in our models range between 0.10 — 
0.11 Mq for Class I objects and 0.06-0.12 Mq for Class 
II objects, well in excess of the AW's median masses for 
either Class I disks (3 x 10~^ -^0) or, especially. Class II 
disks (3 X 10^'^ Mq). Much the same, our Class disks 
have mean mas ses (3.09-0.10 Mq , while five Class disks 
in a sample of Bro wn et al.l (!2000) have a mean mass of 
order 0.01 Mq, though the statistics in the case of Class 
objects is inadequate to draw any firm conclusions. 

A closer look at the data provided by AW reveals that 
our numerical modeling fails to reproduce only the lower 
bound on the measured disk masses but succeeds at re- 
producing the upper bound. For instance, figure 15 of 
AW indicates that Class I disk masses lie in the range 
Md,ci = 0.001 - 0.6 Mq, while our Class I disks have 
masses in the range (Md,ci) = 0.02 — 0.55 Mq. Simi- 
larly, AW's Class II disk masses lie in the range Md,cii = 
5 X 10""* — 0.5 Mq, while our disks have masses in the 
range (A/d,cii) = 4x 10~^— 0.7 Mq. If disk masses are in- 
deed systematically underestimated by conventional ob- 
servational methods, then low-mass disks should suffer 
from this problem to a greater extent than do upper- 
mass disks. The reason for this is not well understood. 

It is worth noting that our Class disk have masses 
that are similar to those of Class I disks. This is in agree- 
ment with the modeling of circumstellar envelo pe dust 
emissi on of 6 deeply embedded systems by Loone v~et al.l 
(|2003f ). who have found that Class systems do not have 
disks that are significantly more massive than those in 
Class I systems. This led them to conclude that the disk 
in Class systems must quickly and efficiently process 
~ 1.0 Mq of material from the envelope onto the pro- 
tostar. Our numerical modeling corroborates their con- 
clusion ~ Class O/I disks develop vigorous gravitational 
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instability that helps keep the disk mass well below that 
of the protostar. 

To what extent can our numerical modeling overes- 
timate disk masses? The reason for the overestimate 
may be twofold. First, Class II disk masses in Figures [2] 
and [5] were derived from disks which have ages of or- 
der 3 Myr but real disks may be older. To test this 
possibility, we ran a few models for 5 Myr and compared 
the resulting Class II disk masses with those obtained for 
3 Myr-old disks. We found that longer ages of disks could 
only cause a factor of 2 (at most) decrease in the time- 
averaged Class II disk masses, still not enough to recon- 
cile our obtained disk masses with those of AW. Second, 
we might have not taken into account some important 
physical mechanisms that help reduce disk masses. In- 
deed, dust disks and accretion disappear on time scales of 
about 6 Myr and this points to the existence of additional 
disk clearing mechanisms, especially relative to the non- 
viscous disk models. Magnetic braking is known to be 
efficient at transporting disk angular momentum directly 
to the external environment. However, disks are known 
to have low ionization fractions and ambipolar diffusion 
may strongly moderate the effect of magnetic braking. 
We plan to explore this mechanism in a follow-up pa- 
per. The formation of a binary system can also decrease 
the resulting disk masses but the fraction of binary sys- 
tems is under debate and it may be too lo w to reconcil e 
our model disk masses with observations (|Ladal [20061 ). 
Photoevaporation of disks due to the external ultravio- 
let radiation may somewhat decrease disk masses but it 
operates only late in the evolution of Class II disks and 
is expected to have little influence on the time-averaged 
disk masses. Clearly, the statistics on Class II disks 
must depend on the existence of additional disk clearing 
mechanism and the adopted end time in our numerical 
and more work is needed to resolve the cause of dispar- 
ity between observationally and numerically derived disk 
masses. 

6. NUMERICAL AND MODEL CAVEATS 

6.1. Gas thermodynamics 

It has recently become evident that the gas thermo- 
dynamics plays an important role in regulating grav- 
itational insta bility and fragment ation in circumstel- 
lar disks (e.g. Pickett et al. 2003; 'Johnson & Gammi^ 
[2003; Rice et al. 2003; Bolcy ct al. 2006; Cai et al. 200^ 
While most researches agree that circumstellar disks are 
susceptible to gravitational instability, particularly in 
the early phase of stellar evolution, there still no con- 
sensus on whether fragmentation ever occurs and, if it 
does, whether any clumps t hat form will become bound 
protoplanets. For instance, Ijohnson fc Gammid ()2003l ) 
have considered geometrically thin disks that cool by ra- 
diatively transporting thermal energy (generated in the 
midplane) to the disk surface. They have found that 
disks are susceptible to fragmentation only if the effec- 
tive cooling time is comparable to or smaller than the 
dynamical time. Simi lar conclusio ns have been made by 
iRice et al.l (|2003f ) and[Meiia| (12004 ) using global SPH and 
grid-based numerical simulations, respectively. Fragmen- 
tation ca n be stabilize d in the inner disks by slow cool- 
ing (e.g. iRafikovl [20051 : iStamatellos fc Whitworthl[200ah 
and in the outer disks by stellar and enve lope irradiation 
(jMatzner fc Levinll200a ICai et all [2001 . On the outer 



hand, disk fragmentation can be aided by the infall of 
matter from an external envelope, certainly in the early 
phase of disk evolution (IVorobvov fc Basul[2005bl 120061 ; 
iKrumholz et al.ll2007l ; fKratter et al.ll2008[ ). 

The accurate implementation of gas thermodynamics 
in numerical simulations of self-consistent disk formation 
and long-term evolution is a very difhcult task. This is 
because our numerical simulations capture several phys- 
ically distinct regimes that see a substantial change with 
time in the chemical composition, opacities, and dust 
properties. Numerical techniques allowing for the accu- 
rate treatment of thermal physics in su ch numerical simu- 
lations are only starting to emerge (e.g. IStamatellos et al.l 
120071)1 ) ■ Taking into account the complexity of gas ther- 
modynamics in circumstellar disks, we use the barotropic 
equation of state. 

It should be stressed here that circumstellar disks de- 
scribed by the barotropic equation of state with 7 = 1.4 
are susceptible to fragmentation and formation of stable 
clumps in the early embedded phase of evolution. These 
clumps are later driven onto the protostar and produce 
bursts of mass accretion that can help explain FU Ori 
eruptions. This phase of protostellar accretion is known 
as the burst phase, it operates in disk-star systems with 
mass > 0.6 Mq and it is very efficie nt at tra nsporting 
the di s k mat erial onto the protostar (IVorobvo v & Bas^ 
I2005bl I2006D . Circumstellar disks described by 7 = 
5/3 are hotter and less susceptible to fragmentation 
but the clump productio n is not completely suppressed 
([Vorobvov fc Basul I2OO6I ). Recent numerical and semi- 
analytic modeling using a more accurate prescription for 
the energy balance in the disk (including radiation trans- 
fer) also predict clump formation in disks (particularly, 
in their outer parts) around stars wit h mass equal to 
or more massive than one solar mass (Krumholz et aLl 
2007; Mayer et al. 2007; StamatcUos et al. 2007a; Bosi 
I2OO8I: iKratter et al. 2008). 

In order to examine if a higher disk temperature can af- 
fect the accretion history and, consequently, disk masses, 
we stiffened the barotropic equation of state by rasing 7 
from 1.4 to 5/3. In the case of self-gravitating disks, 
this affects most of the disk material because most of 
the disk is optically thick and is characterized by surface 
densities considerably la rger than Scrit = 36.2 g cm^^ 
(| Vorobvov fc Basul [20071) . On the other hand. Class II 
viscous disks may have surface densities comparable to or 
lower than Scriti pa rticularly in the late evolution phase 
([Vorobvov fc Basul [2008b ). Even in this case, the in- 
crease in 7 is expected to affect a large portion of the disk 
material because equation ([3]) makes a gradual transition 
between the optically thin E < < Scrit and optically thick 
E >> Ecrit regi mes. The detailed numeric al simulations 
are presented in IVorobvov fc Basul (|2008b[ ). here we pro- 
vide only the main results. We find that while an increase 
in 7 raises the disk temperature from 50 K (at 10 AU) to 
about 100 K and makes the disk less susceptible to frag- 
mentation and clump formation, the amount of accreted 
material changes insignificantly. In fact, the mass accre- 
tion rates time-averaged over ~ 10* yr are very similar 
in models with 7 = 1.4 and 7 = 5/3, while the instan- 
taneous rates may be quite different. This is because 
a modest increase in disk temperature acts to suppress 
higher order spiral modes m > 2, which are rather in- 
efficient at transporting mass radially inward and tend 
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to produce more fluctuations and cancellation in the net 
gravitational torque. At the same time, the growth of 
lower order modes m < 2 is promoted, which are efficient 
agents for mass and angular m omentum r edistr ibution. 
A similar effect was reported by ICai et aH (1200 8*) for cir- 
cumstellar disks heated by a strong envelope irradiation. 
Thus, the relative strength of lower order spiral modes 
increases in hotter disks, which compensates an appar- 
ent decrease in the amount of accreted gas due to the 
less efficient burst phase of accretion. 

6.2. Numerical resolution and central point mass 

The use of a logarithmically spaced numerical grid in 
the radial direction means that the ratio Ar/r of the cell 
size Ar to radius r is constant for a given cloud core and 
varies from 0.05 (model 21) to 0.07 (model 7). It can 
be noticed from Fig. lAll that the aspect ratio A ^ z/r 
of the disk vertical scale height to radius is smaller than 
Ar/r in the inner 70 AU, which implies that a better 
radial (and angular) resolution is desirable in this inner 
region. We have found that an increase in the resolu- 
tion to 256 X 256 grid points makes little influence on the 
accretion history and disk masses but helps save a con- 
siderable amount of CPU time and, eventually, consider 
many more model cloud cores. We also point out that, 
due to the adopted log spacing in the radial direction, the 
resolution in the inner computational region is sufficient 
to fulfill the Truelove criterion in the disk. 

In our numerical simulations the position of the central 
star is fixed in the coordinate center. In reality, however, 
the star may wobble around the center of mass in re- 
sponse to the gravity force of the disk. This can promote 
the growth of non-axisymmetric spiral modes (especially 
the m=l mode) in massive disks and, consequently, help 
limit disk masses in the early phas e of stellar evolution 
(|Adams et al.lll989l : IShu et ai.lll990f ). We plan to explore 
this potentially important mechanism in a follow-up pa- 
per. 

6.3. The weight function 

The mean disk masses listed in Table |4] may depend 
on the form of the weight function. Our adopted weight 
function J-'k [see eq. (jlip ] assumes that stellar masses in 
each stellar evoluti on phase are d istributed according to 
the Kroupa IMF ([Kroupal [2001h . which actually refers 
to the initial masses of stars that have already formed. 
To determine the extent to which our mean disk masses 
depend on the particular form of the weight function, we 
also consider weight functions with a slop e typical for th e 
Salpeter IMF, Ts ((M.. )) = A(M.)-^-^-^ (iSalpete d flQSl . 
and MiUer-Scalo IMF (jMiller fc Scalolfl979[ ) 

( A (Af,)-i-25 if 0.1 Mq < (M,) < 1.0 M( 
^MS ((M*)) = < B (Af,)-2-0 if 1.0 Mq < {A-Q < 2.0 M( 
( C (Af,)-2-3 if (^M,) > 2.0 Mq. 

(17) 

The resulted mean disk masses are listed in Table [3 It 
is evident that the use of different weight functions (in- 
dicated in parentheses) yields mean values that are dif- 
ferent from those derived from the Kroupa law by 30% 
at most. 

7. SUMMARY 



TABLE 5 

Mean disk masses with different weight functions 



disk type 


A^d.CO 


A/d,ci 


A^d.CII 


self-gravitating (J^s) 


0.08 


0.08 


0.10 


viscous (J^s) 


0.09 


0.09 


0.05 


self-gravitating (J-'ms) 


0.09 


0.12 


0.15 


viscous (J^Ms) 


0.10 


0.13 


0.08 



Note. — All disk masses are in Mq. The weight functions used 
to derive the mean disk masses are shown in parentheses. 

We have considered numerically the long-term evo- 
lution of self-consistently formed self-gravitating and 
viscous circumstellar disks around low-mass stars 
(0.2 Mq < < 2.0 Mq). We seek to determine 
time-averaged disk masses in the Class ((Md^co)), 
Class I ((Md,ci)), and Class II ((Md,cii)) stellar evo- 
lution phases. Time averaging in each evolution phase 
is done over the duration of the phase. The mean disk 
masses (Md) are then derived from these time- averaged 
masses by weighting them over the corresponding stellar 
masses using a power-law function with a slope typical 
for the Kroupa initial mass function of star (see eq. [12]). 
In self-gravitating disks the radial transport of mass and 
angular momentum is performed exclusively by gravita- 
tional torques, while in viscous disks both the gravita- 
tional and viscous torques are at work. Our numerical 
simulations yield the following results. 

1. Both the self-gravitating and viscous disks have 
Class I mean masses Md,ci = 0.10 — 0.11 Mq that 
are slightly more massive than those of Class 
disks Md,co = 0.09 - 0.10 Mq. However, vis- 
cous Class II disks have a mean mass (Md,cii = 
0.06 Mq) that is a factor of 2 lower than that of 
self-gravitating Class II disks. This is explained by 
the fact that gravitational torques prevail through 
the Class and Class I phases but succumb to vis- 
cous torques through much of the Class II phase. 

2. Our obtained mean disk masses are larger than 
those derived by AW for 153 YSOs in the Taurus- 
Auriga star formation region, regardless of the 
physical mechanisms of mass transport in the disk. 
The difference is especially large for Class II disks, 
for which we find Md,cii = 0.06-0.12 Mq but AW 
report median masses of order 3 x 10"'^ Mq. Our 
Class I disks have almost a factor of 4 larger mean 
disk mass than those of AW. 

3. Time-averaged disk masses have a considerable 
scatter around mean values in each evolution phase. 
For instance. Class and Class I disk masses 
lie in the range (Md,co) = 0.04 - 0.2 Mq and 
(A^d,ci) = 0.02 - 0.55 Mq, respectively, while 
Class II disks have a much wider range in masses 
Md,cii = 0.004-0.7 Mq. 

4. When the three stellar evolution phases are con- 
sidered altogether, we find a near-linear relation 
between time-averaged disk and stellar masses, 
(Md) = (0.2 ± 0.05) (Mh.)1-3±o i5. This rela- 
tion can potentially become somewhat shallower if 
cloud cores with very low ratios of rotational-to- 
gravitational acceleration are considered, but can- 
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not be broken completely. However, when each 
phase is considered separately, the relation be- 
tween disk and stellar masses becomes progres- 
sively steeper along the sequence of stellar evolu- 
tion phases. In particular, the corresponding re- 
lations for Class 0, Class I, and Class II objects 
have exponents 0.7 ± 0.2, 1.3 ± 0.15, and 2.2 ± 0.2, 
respectively. 

In each stellar evolution phase, time-averaged disk- 
to-star mass ratios {£,) tend to have greater val- 
ues for stars of greater mass. However, there is a 
clear saturation effect - (^) never exceed 40%, re- 
gardless of the stellar mass. This means that no 
star can harbour a disk with mass comparable to 
or greater than that of the star. This saturation ef- 
fect is caused by the onset of vigorous gravitational 
instability in circumstellar disks tha t form in the 
Class or early Class I phase (see also lAdams et al.l 
[l98l[siurerSlT990l) . 



Only low-mass objects < 0.9 Mq are expected 
to have Class disks. Some Class objects (par- 
ticularly those formed from slowly rotating cloud 
cores) may have no disks and outflows associated 
with them. 
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APPENDIX 
DISK SCALE HEIGHT 

We derive the disk vertical scale height Z at each computational cell via the equation of local vertical pressure 
balance 



pel 



■ 9 Z.St 



(Al) 



where p is the gas volume density, gz,gas and g^^st are the vertical gravitational accelerations due to self-gravity of a 
gas layer and gravitational pull of a central star, respectively. Assuming that p is a slowly varying function of vertical 



distance z between z 



(midplane) and z 

P9z 



Z (i.e. S = 2 Z p) and using the Gauss theorem, one can show that 



i dz : 



pgz,st dz^ 
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where r is the radial distance and Af, is the mass of the central star. Substituting equations (jA2p and (jA3P back into 
equation (jAip we obtain 
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2GAUp 
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(A4) 



This can be solved for p given the model's known c^, S, and using Newton-Raphson iterations. The vertical scale 
height is finally derived as Z = I]/(2p). 

Finally, we compare our model values o f Z with those p r edicte d from detailed vertical structure models of irradiated 
accretion disks around T Tauri stars bv lD'Alessio et al] (119991 ). Their figure lb (dashed curve) yields the following 
relation between the aspect ratio A = Z/r of the local scale height to radial distance and radial distance r 

A = 0.1(r/100AU)^/\ (A5) 
which is plotted by the dashed line in Figure [XT] The solid line in Figure ET] shows our obtained aspect ratio A = Z/r 
for model 13 at t=0.3 Myr. It is evident that our A is somewhat underestimated in the inner disk and overestimated 
in the outer disk but the disagreement between the two curves is within a factor of unity. 



DIVERGENCE OF THE VISCOUS STRESS TENSOR 



The components of V • 11 in polar coordinates (r, 0) are 

1 d 

(V-n)^ = -— rH,, 



dr 

(v-n), = ^H 



ld_ 

r d(f> 



H, 



+ 2- 



dr""^'' ^ r dct^"^^ ~^ " r 
where we have neglected the contribution from off-diagonal co mponents Hrz and H, 



-rm^^ 



Hr 



(Bl) 
(B2) 



stress tensor 11 in polar coordinates (r, i 



The components of the viscous 

can be found in e.g. IVorobvov fc Theii ()2006l ). 
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